In [8]:
def train_and_test_model(df, xvars, yvar,
to_plot=1, model=None, num_experiments=None,
log=False, square=False, sqrt=False):
"""
Train and test a univariate linear regression model with cross-validation
Args:
df is a Pandas DataFrame object
xvars is a list of strings (the feature names)
yvar is a string (the label name)
to_plot : make debug plots if true
Returns:
a dictionary holding model outputs
"""
# Fraction of samples reserved for model testing
# e.g. 0.2 means 20% of samples reserved. In the
# context of cross-validation, this would be
# equivalent to a 5-fold cross-validation. The
# cross-validation we end up doing is not exactly
# 5-fold cross-validation since we introduce random
# sampling and many many repeats
TEST_SIZE = 0.2
# Number of trials
num_experiments = num_experiments or 100
# Unpack and filter data
xvars = [xvars] if isinstance(xvars, str) else xvars
# Create a temporary df to filter NaNs only from the columns we use
xyvars = xvars + [yvar]
df_nonans = df[xyvars].dropna()
if log is True:
df_nonans = df_nonans.apply(lambda x: np.log(abs(x)))
if square is True:
df_nonans = df_nonans.apply(lambda x: x**2)
if sqrt is True:
df_nonans = df_nonans.apply(lambda x: np.sqrt(abs(x)))
X = df_nonans[xvars].values
y = df_nonans[yvar].values
percent_error_list = []
percent_error_train_list = []
r2_score_list = []
rmse_list = []
baseline_rmse_list = []
intercept_list = []
coef_list = []
corr_list = []
has_plotted = 0
for idx in range(num_experiments):
X_train, X_test, y_train, y_test = train_test_split(X, y,
test_size=TEST_SIZE)
# apply standard scaler (fit only on train, transform both train and test)
# keep originals for pltos
X_train_original = X_train
X_test_original = X_test
scaler = StandardScaler()
X_train = scaler.fit_transform(X_train)
X_test = scaler.transform(X_test)
(model, corr) = build_model(X_train, y_train, model)
y_pred_train = model.predict(X_train)
y_pred = model.predict(X_test)
y_train_mean = np.mean(y_train)
r2_score = model.score(X_test, y_test)
if log is True:
y_train = np.exp(y_train)
y_pred_train = np.exp(y_pred_train)
y_pred = np.exp(y_pred)
y_test = np.exp(y_test)
y_train_mean = np.exp(y_train_mean)
if square is True:
y_train = np.sqrt(y_train)
y_pred_train = np.sqrt(y_pred_train)
y_pred = np.sqrt(y_pred)
y_test = np.sqrt(y_test)
y_train_mean = np.sqrt(y_train_mean)
if sqrt is True:
y_train = (y_train)**2
y_pred_train = (y_pred_train)**2
y_pred = (y_pred)**2
y_test = (y_test)**2
y_train_mean = (y_train_mean)**2
rmse = np.sqrt(np.mean((y_pred - y_test)**2))
percent_error = np.mean(np.abs(y_pred - y_test)/y_test) * 100
baseline_rmse = np.sqrt(np.mean((y_train_mean - y_test)**2))
percent_error_train = np.mean(np.abs(y_pred_train - y_train)/y_train) * 100
rmse_list.append(rmse)
percent_error_list.append(percent_error)
percent_error_train_list.append(percent_error_train)
r2_score_list.append(r2_score)
baseline_rmse_list.append(baseline_rmse)
# intercept_list.append(model.intercept_)
# coef_list.append(model.coef_)
corr_list.append(corr)
# Give a visual into the cross-validation process
if to_plot and not has_plotted and X.shape[1] == 1:# and baseline_rmse < rmse:
# print(r2_score, rmse, baseline_rmse)
has_plotted = 1
x_fit = np.linspace(np.min(X), np.max(X), 100)
x_fit_trans = scaler.transform(x_fit.reshape(-1,1))
y_fit = model.predict(x_fit_trans)
if log is True:
y_fit = np.exp(y_fit)
if square is True:
y_fit = np.sqrt(y_fit)
if sqrt is True:
y_fit = y_fit ** 2
plt.plot(x_fit, y_fit, color=(0.3, 0.3, 0.3))
plt.plot(X_train_original, y_train, 'o', color='k')
plt.plot(X_test_original, y_test, 's', color='r')
for (this_x, this_y, this_ypred) in zip(X_test_original, y_test, y_pred):
plt.plot([this_x, this_x],
[this_y, this_ypred],
label='', color=(0.3, 0.3, 0.3))
plt.xlabel('x')
plt.ylabel('y')
plt.legend(['Model', 'Train', 'Test'])
plt.show()
# Package the output
final_model = model
# Column mean and column std of a list of tuples
# final_model.coef_ = np.array(list(map(np.mean, zip(*coef_list))))
# final_model.intercept_ = np.mean(intercept_list)
final_rmse = np.mean(rmse_list)
final_rmse_std = np.std(rmse_list)
final_baseline_rmse = np.mean(baseline_rmse_list)
final_baseline_rmse_std = np.std(baseline_rmse_list)
# Column mean and column std of a list of tuples
final_corr = tuple(map(np.mean, zip(*corr_list)))
final_corr_std = tuple(map(np.std, zip(*corr_list)))
final_percent_error = np.mean(percent_error_list)
final_percent_error_std = np.std(percent_error_list)
final_percent_error_train = np.mean(percent_error_train_list)
final_percent_error_train_std = np.std(percent_error_train_list)
# print(len(r2_score_list), sum(np.array(r2_score_list)>0))
final_r2_score = np.mean(r2_score_list)
final_r2_score_std = np.std(r2_score_list)
model_output_dict = dict()
model_output_dict['model'] = final_model
model_output_dict['rmse'] = final_rmse
model_output_dict['rmse_std'] = final_rmse_std
model_output_dict['baseline_rmse'] = final_baseline_rmse
model_output_dict['baseline_rmse_std'] = final_baseline_rmse_std
model_output_dict['corr'] = final_corr
model_output_dict['corr_std'] = final_corr_std
model_output_dict['percent_error'] = final_percent_error
model_output_dict['percent_error_std'] = final_percent_error_std
model_output_dict['percent_error_train'] = final_percent_error_train
model_output_dict['percent_error_train_std'] = final_percent_error_train_std
model_output_dict['r2_score'] = final_r2_score
model_output_dict['r2_score_std'] = final_r2_score_std
model_output_dict['r2_score_manual'] = 1 - final_rmse ** 2 / final_baseline_rmse ** 2
model_output_dict['features'] = tuple(xvars)
return model_output_dict
# Test
# one feature
# train_and_test_model(df,
# "var_q_c20_c365_c3_ah",
# 'cycles_to_80_pct',
# to_plot=1)
# train_and_test_model(df,
# "var_q_c20_c365_c3_ah",
# 'cycles_to_80_pct',
# to_plot=1,
# log=True)
# train_and_test_model(df,
# "var_q_c20_c365_c3_ah",
# 'cycles_to_80_pct',
# to_plot=1,
# square=True)
train_and_test_model(df,
"var_q_c20_c365_c3_ah",
'cycles_to_80_pct',
to_plot=1,
sqrt=True)
# multiple features
train_and_test_model(df,
["var_q_c20_c365_c3_ah", "var_q_c20_c56_c3_ah", "var_q_c20_c159_c3_ah"],
'cycles_to_80_pct',
to_plot=1,
model=ElasticNet())
Out[8]:
{'model': ElasticNet(),
'rmse': 25.02531852609627,
'rmse_std': 12.317017918591626,
'baseline_rmse': 42.54946122739829,
'baseline_rmse_std': 13.694390349187083,
'corr': (-0.8769411874917533, -0.5812288664987114, -0.5304186200087936),
'corr_std': (0.028429759667883758, 0.08068569908181904, 0.12185106756653163),
'percent_error': 5.244552038299888,
'percent_error_std': 2.1709569835195506,
'percent_error_train': 4.4514011570488945,
'percent_error_train_std': 0.4656297602449225,
'r2_score': 0.18524711379733433,
'r2_score_std': 1.2024225388716987,
'r2_score_manual': 0.6540836439933803,
'features': ('var_q_c20_c365_c3_ah',
'var_q_c20_c56_c3_ah',
'var_q_c20_c159_c3_ah')}
